1 Poisson Events

1.1 Libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))

1.2 Parameters and risk

censoredProb <- 0.0002
timeSpan <- 20
timeInterval = 0.01
InitialPopulatoin <- 1000
ContBetaRate_1 <- 0.0002
ContBetaRate_2 <- 0.000001
BinBetaRate_1 <- 0.001
BinBetaRate_2 <- 0.003
BaselineHazard <- ContBetaRate_1
betaRates <- c(BaselineHazard,ContBetaRate_1,ContBetaRate_2,BinBetaRate_1,BinBetaRate_2)
BaseVar <- rep(1,InitialPopulatoin)
ContVar_1 <- runif(InitialPopulatoin)
summary(ContVar_1)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001145 0.2347710 0.4894608 0.4923210 0.7424088 0.9995301
ContVar_2 <- rnorm(InitialPopulatoin,50,10)
ContVar_2[ContVar_2 < 1] <- 1 
summary(ContVar_2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.41   42.91   49.26   49.62   55.89   83.67
BinVar_1 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_1)
## BinVar_1
##   0   1 
## 639 361
BinVar_2 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_2)
## BinVar_2
##   0   1 
## 648 352
dataFeatures <- as.matrix(cbind(BaseVar,ContVar_1,ContVar_2,BinVar_1,BinVar_2))
hazardRate <- as.numeric(dataFeatures %*% betaRates)
summary(hazardRate)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002355 0.0003632 0.0013097 0.0017651 0.0033525 0.0044680

1.2.1 Getting the events and time to event

aliveSet <- c(1:InitialPopulatoin)
eventSet <- numeric(InitialPopulatoin)
timetoEvent <- numeric(InitialPopulatoin)

for (time in c(1:(timeSpan/timeInterval)))
{
  randProb <- runif(length(aliveSet))
  Iscensored <- randProb <= censoredProb
  Isevent <- randProb <= (1.0-exp(-hazardRate[aliveSet]))
  Iscensored <- Iscensored & !Isevent
  eventSet[aliveSet] <- Isevent
  timetoEvent[aliveSet] <- time*timeInterval-timeInterval/2
  isCensoredOrEvent <- Iscensored | Isevent
  aliveSet <- aliveSet[!isCensoredOrEvent]
#  cat(length(aliveSet),"(",sum(isCensoredOrEvent),",",sum(Isevent),",",sum(Iscensored),")\n")
}
timetoEvent[aliveSet] <- time*timeInterval + timeInterval/2
summary(timetoEvent)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.005   1.603   5.040   8.283  16.183  20.005
table(eventSet)
## eventSet
##   0   1 
## 203 797
pevent <- (1.0-exp(-hazardRate))
summary(pevent)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002355 0.0003631 0.0013088 0.0017624 0.0033469 0.0044580
simulatedDataFrame <- as.data.frame(cbind(status=eventSet,time=timetoEvent,pevent=pevent,dataFeatures))

1.3 RRplots()

plotTimeInterval <- 2.0

hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)


rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04600 0.07007 0.23044 0.26640 0.48855 0.59082

table(simulatedDataFrame$status)

0 1 203 797

RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=simulatedDataFrame$time,
                     title="Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=plotTimeInterval)

par(op)

1.3.1 By Risk Categories

obsexp <- RRAnalysisCI$OERatio$atThrEstimates
pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 797 742.6 854.3 758.2 1.05 0.979 1.13 0.162
low 190 163.9 219.0 162.4 1.17 1.010 1.35 0.034
90% 32 21.9 45.2 29.3 1.09 0.746 1.54 0.579
80% 575 529.0 624.0 567.9 1.01 0.931 1.10 0.753 
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))

plot(obsexp$Expected,obsexp$Observed,
     xlim=c(minx,maxx),
     ylim=c(minx,maxx),
     main="Cal. Expected vs Observed",
     ylab="Observed",
     xlab="Expected",
     col=rainbow(nrow(obsexp)),
     log="xy")

errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

1.3.2 Time to Event Analysis

isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Observed Time",
                   at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames") 
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
  O:Low Risk < 0.080 O:0.080 <= Risk < 0.085
1Q 3.26 3.89
Median 8.03 8.43
3Q 14.18 14.26
Table continues below
  O:High Risk >= 0.085 E:Low Risk < 0.080
1Q 1.02 13.3
Median 2.35 14.8
3Q 5.34 17.0
  E:0.080 <= Risk < 0.085 E:High Risk >= 0.085
1Q 11.4 1.45
Median 11.6 1.51
3Q 11.8 3.65 
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
     xlab="Mean Expected Time",
     ylab="Mean Observed",
     xlim=c(minv,maxv),
     ylim=c(minv,maxv),
     main="Estimated Time to Event",col=rainbow(length(classnames)))

lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

1.3.3 Risk Calibration

op <- par(no.readonly = TRUE)


crdata <- cbind(simulatedDataFrame$status,pvalue,simulatedDataFrame$time)

#calprob <- CalibrationProbPoissonRisk(crdata,timeInterval=plotTimeInterval)
calprob <- CalibrationProbPoissonRisk(crdata)

( 20.72695 , 11020.74 , 1.956533 , 797 , 1588.147 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
7.24 11.1 20.7

1.3.4 After Calibration

h0 <- calprob$h0

caldata <- cbind(simulatedDataFrame$status,calprob$prob)


rrAnalysisTrain <- RRPlot(caldata,atRate=c(0.90,0.80),
                     timetoEvent=simulatedDataFrame$time,
                     title="Cal. Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

1.3.5 By Risk Categories

obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 797 742.6 854.3 811.7 0.982 0.915 1.05 0.623
low 190 163.9 219.0 173.8 1.093 0.943 1.26 0.225
90% 32 21.9 45.2 31.4 1.019 0.697 1.44 0.858
80% 575 529.0 624.0 607.9 0.946 0.870 1.03 0.187 
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))

plot(obsexp$Expected,obsexp$Observed,
     xlim=c(minx,maxx),
     ylim=c(minx,maxx),
     main="Cal. Expected vs Observed",
     ylab="Observed",
     xlab="Expected",
     col=rainbow(nrow(obsexp)),
     log="xy")

errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

1.3.6 Time to Event Analysis

isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Observed Time",
                   at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames") 
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
  O:Low Risk < 0.080 O:0.080 <= Risk < 0.085
1Q 3.26 3.89
Median 8.03 8.43
3Q 14.18 14.26
Table continues below
  O:High Risk >= 0.085 E:Low Risk < 0.080
1Q 1.02 13.3
Median 2.35 14.8
3Q 5.34 17.0
  E:0.080 <= Risk < 0.085 E:High Risk >= 0.085
1Q 11.4 1.45
Median 11.6 1.51
3Q 11.8 3.65 
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
     xlab="Mean Expected Time",
     ylab="Mean Observed",
     xlim=c(minv,maxv),
     ylim=c(minv,maxv),
     main="Estimated Time to Event",col=rainbow(length(classnames)))

lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)